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I. INTRODUCTION: 


1.1 General Statement of the Problem 

On the event of loss of vacuum guard of superinsulated helium dewar high 
rate of heat transfer into the tank occurs. Rapid boiling of liquid helium 
causes burst disk to rupture at four atmospheres and consequent passage of 
helium to atmosphere through the vent lines. Gaseous helium exiting the vent 
line forms vertical buoyant jet in a stagnant environment. 


Characterization of the gaseous jet is achieved by detailed analysis of 

I 

axial and radial dependence of the flow parameters. Unsteady flow pattern at 
the jet exit influences the developing profile downstream. Figure 1.1 is a 
schematic representation of the asymmetric jet. Three identifiable regimes of 
the jet are illustrated in the figure. Such trend of the radial profiles are 
observed in constant density jets and can also be adapted for variable density 
jet through some corrections in effective jet parameters. The potential core 


is a part of the developing shear layer where the jet is assumed to be 


uniform 


and inviscidy jf jtelocity^ temperature, density and concentration at the axis of 
the potential core is assumed^ to be equal to jet discharge values. Pressure 
recovery from the upstream throat pressure to atmospheric pressure occurs in 
this region. / 


The process of jet development include several complex phenomenon 
including turbulence, while overall character of the jet is determined by the 
strength of the global forces effective in the fluid motion. The fluid motion 
in a buoyant jet is in general governed by buoyant, viscous and inertial 
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forces. The relative magnitude of these forces define the local character of 
the jet in different regions. On the other hand, the overall character is 
determined by the magnitude of these forces existing at the jet source and the 
ambient condition. Exit Reynolds number, Froude number and Grashof number are 
generally used to characterize the jet as explained in section 2.1. 

The final regime illustrated in figure 1.1 is called the "self-similarity 
regime./ A flow field can be called self-similar. when only one geometrical 
variable is required to characterize the nondimensional ized, time averaged 
behavior-of velocity, temperature or concentration throughout the region. 
Complete self-simi larityiY never achieved by a developing flow. Yet, in an: 
analogous way, local self-similarity can be defined through local geometric 
variables.^ For example, the radius of U max /2 can play the role of local 
length scale used to normalize the radial distribution functions. 

Entrainment of surrounding air into the jet causes broadening of the jet 
width and creates a mixture zone. Since vaporization of helium is completed 
upstream from the jet nozzle, single phase analysis is adequate to describe 
the dynamics of the jet. On the other hand, excessive variation of tem- 
perature induces predominant variation of gas density. Hence variable density 
single phase model has been considered here. 

Theoretical considerations given to helium jet dispersion analysis are 
illustrated in Chapter II. The results and analysis obtained from the com- 
puter program developed for this project are illustrated in section 3.1. 
Prediction of axial and radial distribution of temperature and velocity are 
emphasized for illustraton. Axial velocity distributions indicate nonlinear 


2 


decay profiles while axial temperature distributions indicate asymptotic 
increase from jet exit to surrounding temperature. Additional analysis using 
a second solution scheme is performed with computer program GENMIX. The 
results of this analysis are illustrated in section 3.2. Single phase helium 

a 

* 

jet analysis is considered with identical initial and boundary conditions 
used in both analysis. Jet discharge conditions are variable parameters to 
the programs. 

An user's guide for the HEJET program is provided in appendix A. Input 
and output variable sequences are illustrated using corresponding file descrip- 
tions. , 
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1.2 Previous Related Studies 


Jet dispersion analysis has been studied widely under different hydrodyna- 
mic and thermodynamic circumstances. Research interest has focused on single 

* 

and two phase jet characterization through effects of temperature, velocity and 
jet to ambient density ratios. A detailed review of the experimental data on 
single vertical buoyant jets are given by Chen and Rodi (1). Similarity and 
scaling laws are discussed by the author for jet characterization. Axial 
and radial distribution of mean velocity, temperature and concentration of 
jets are compared for a wide range experiments. 

t 

Two phase jet measurements with emphasis in spray evaporation has been 
reviewed by Shearer and Faeth (2). The authors have considered evaporation of 
a well-atomized liquid jet where homogeneous equilibrium model has been 
assumed and validated. Second order turbulence modeling including fluctuating 
density effects was employed for detailed analysis. The authors have also 
produced measurements of single phase variable density jets using sulfur 
hexafluoride gas jet. Freon-11 spray was produced by air atomizing injector 
in order to produce an evaporative jet. 

A complete review of jet characterization measurements is given by Pitts 
(3). Emphasizing the effect of the ratio of jet fluid density to surrounding 
air density, the entire spectrum of experiments were suiunarized. Velocity 
of the similarity analysis was demonstrated for majority of the experiments. 
Rayleigh scattering technique was utilized for flow visualization and cen- 
terline concentration measurements. Gases of wide range of density were con- 
sidered for analysis and compared to existing data. 
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Helium jet analysis is shown in figure 1.2. Characteristics of rate of 
spread of the jet has been demonstrated through the similarity parameters. 


K 
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II THEORETICAL CONSIDERATIONS . 




ft 


> 
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2.1 Description of the model 


Similarity analysis has been employed in this study as a means of solving 
the characteristic flow parameters of a buoyant helium jet. With the usual 
boundary layer approximations, the differential equations governing the mean 
flow quantities in a vertical buoyant jet can be written as 


Continuity Equation: 


3(pUr) + 9(pVr) = Q 
9Z 8r 


Momentum Equation: 


+ . . g(p . Pa)r -f^, 


Thermal Energy Equation: 




Concentration: 


3(pUCr) . 3(pVCr) 
3Z + 3r 


= -^(rpvc 1 ) 


( 1 ) 


(2) 


(3) 


(4) 


Integration of equations 1 to 4 over the jet cross-section area 
yields the integral form these equations: 

Continuity Equation ^ 


E 

p Ur dr = E (5) 
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fT , 


Momentum Equation 


d_ 

dZ 


•p U 2 r dr = g 


(p a -p)r dr 


( 6 ) 


r r r 


Thermal Energy 


d_ 

dZ 


dT a 

P U(T-T a )rdr = - ^ 


pUrdr 


(7) 


r r r 


Concentration: 


d 

dZ 


pU(C-C )rdr 

a 


dC. 

C 

dZ 


rr r 


pllrdr 


( 8 ) 


Here E refers to entrained mass per unit length of the jet (divided by 2 (tt)) 
and r £ is the radius of the'jet from the axis of symmetry. 

In order. to solve the set of equations (5) (6) (7) & (8), certain 
non-dimensional ized local similarity parameter is specified, such that ; 

^n ; = :;r/(rJ$) u (9) 

where . (rJ4) u ; = is the radius of U max /2.::. Using the local 
similarity parameter n, it can be shown. 


c> c m fj (0) (10) 

and 

u = u m f 2 <n) ; ( 11 ) 

where C and U refer to the maximum values of mean concentration and 
m m 

velocity in the radial profiles. Assuming Gaussian distribution. 


(n) = e- -r 


.693 


and 


.693q 2 

U (n) = e" - J 


where R = ^ r ^c 

:: WU 

(riil = radius of C m /2 

_ c - i — . m - 


radius T/2 

m • 
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Implementing equations (10) and (11) into equations (5) - (8), and noting 

HP HT 

that a = 0 and _a = 0 in non-stratified environment for gaseous jets, 
dZ dZ 

we have general ized jet fluid concentration 


C C a 
C o - C a 


T T a 
T o" T a 


p o p 


p p o ‘ p a 


(a) Integral momentum equation 


I 4 _ 


( ~ 7 


1 f 0° ~ f o f 


/ 

'o 


( 12 ) 


M, 


2trp a (r*)J 


U m 

m 


f 2 2 n 


dn 


n 1 + C JPjP n - 1) 
o m' r a 'o ' 1 


(13) 


irr* p n U* 
o r o o 


= momentum of the jet at discharge J 


(b) Integral mass deficit equation 


N. = 2tt p a (rtf)?, U C 
o ' 'u m m 


= irr* C U 
ooo 


fj f 2 n dn 


1 + c m ‘Pa'Po - ') f l 


= Mass deficit of the jet at discharge 


(c) Entrainment law 


(14) 


Vi 

m = k Z M p n = 2™ (r^)„ U m 

e e 0^0 K a ' 'u m 


f 2 n dn 

1 + W»o - » f l 


(15) 


The left hand side of equation (15) has been obtained from measurements 
of Referenced, and confirmed by dimensional analysis. 

Equations (13), (14), and (15) are hence to be solved simultaneously in/ 

order to obtain the unknown parameters U , C and (r&)„ which define the 

. - m m Zju - 


0^ 
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characteristics of the axial profiles. The radial distribution of the flow 
parameters are then achieved through equations (10), (11) and (12) respec- 
tively. 
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2.2 Thermodynamic Properties of Helium 


Complete summary of useful thermodynamic and physical properties of helium 
are available in Reference 5. Considering that the gaseous jet at tains ..... --J 

i ... ....... " 7 . 

atmospheric pressure immediately after the initial shock progression, the ana-> 


lysis has been performed at f atmosphenfc presslTre. Only the discharge tern 


perature therefore influences the. physical and thermodynamic characterics of 
the helium jet. 


Linearity of equation of state is well documented for wide range of 
pressures (5). 

Table 2.1 Thermodynamic and physical properties of helium 


Helium Jet 

Discharge 

Property 

Case 1 

Case 2 

Temperature 

12 K 

120'IC 

Density 

4.0523 ^ 
m J 

0.40523 ^ 
m 3 

Pressure 

1 bar 

1 bar 

Viscosity 

2.23x10-* 

poise 

114x10-* 

poise 

p 

critical 

2.3 bar 

2.3 bar 

^critical 

Mass flow 
rate 

5.2 K 
°.8 

sec 

5.2 K 
°- 8 

sec 





2.3 HEJET Computer Program Description 


> 




H 


The HEJET computer program is written to solve equations (13), (14) and 
(15) simultaneously in order to obtain the axial distributions of centerline 
velocity (U m ) , concentration (C m , and hence temperature T m and density) in 
addition to radius of U m /2 and T m /2 . Subsequently, the program computes 
radial distributions of concentration (temperature or density) and velocity 
from equations (10) and (11). 


The- structure of the program is shown in the tree diagram of Table A.l. 

Function of input and output files are illustrated in Appendix A. Here, the 

« 

main program and subroutine AINTG1 will be discussed. 


Calculation begins with evaluation of the dimensionless numbers such as 
Reynolds number, Froude number and Grashof number. Implication of the 
magnitudes of these numbers are discussed here. The dimensionless forms are: 

(a) Reynolds number R g = U °P 0 D ^- Inertia force ^ 

7 . Viscous force 

^o 


(b) Froude number F 


(c) Grashof number G p 


U 2 p lnertia force -. 

~n7~ — - . 7— “-Buoyant force- 1 

9D(P a -p 0 ) / P 0 


g(p -p )D* r Buoyant force -. 
| — “-Viscous force- 1 
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Purely 


turbulent jets jhave high Reynolds number caused by increased 


inertia relative to viscous effects. On the other hand, pure (turbulent plumes 


are characterized by high Grashof number due to buoyancy e ffects. a( turbulent" 
buoyant .jet] is a combined effect of the two which is characterized by high 
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(^ jnd higl^jTTpk suiting in intermediate value of Froude number (0 < F < ■*»). 

This type of jets are creat ed by discharging fluid of density\7owe?^ than the 
[^density environment./ 


The program then calculates thelpotential core [which is defined as the 


region dominated by inertia and characterized by uniform jet velocity equal to 


V 


The length of the potential cone is estimated by 


x c = 2.13 D (R e ) 


*”1 • 


If the integrals of equation (13), (14) and (15) are denoted by I ^ , I- and 
I 3 respectively, then simultaneous solution of these equations provide 





(16) 

(17) 

(18) 


Subroutine AINTG1 performs the integrals Ij, and using the 

supplementary subroutines QROMO, POLINT, FUNC1, FUNC2, and FUNC3. . Here QROMO 

performs Romberg integration on an open interval where one of the limits may 

extend to infinity. POLINT is used to obtain polynomial of the functions 

FUNC1, FUNC2, and FUNC3 representing the functional forms of I ^ , I 2 and I 3 

respectively. The iterative method is continued until convergence 

i+1 i i 

of C is reached, such that C - C < e C , where e is the desired con- 
m mm m 

vergence criteria. 
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Ill RESULTS AND ANALYSIS 

Self similarity of turbulent buoyant jet has been assumed in the present 
analysis using HEJET computer code. The theoretical background and solution 
scheme has already been discussed in Chapter II. Self similarity of the 
profiles imply that the velocity U and the temperature T can be expressed as 


where n = r/(r&)u * 

subscript ni' "indicates maximum radial value ap pearing at jet axis and subscript 
^corresponds to atmospheric condition. 

Accordingly, it is assumed that the dimensionless form of time averaged 
quantities of a two-dimensional (ax i symmetric) buoyant jet can be well 
described by a single normalized length parameter n- By incorporating the 
above distributions of velocity and temperature into mass, momemtum and energy 
balance equations (1,2 and 3 respectively), simultaneous solution of axial 
distribution of U^, T m and (r!4) u are performed. The asymptotic value of jet 
temperature approaches the ambient temperature, which is assumed to be 
stagnant. Radial distributions of velocity and temperature are thereby calcu- 
lated using the distribution functions fj, and f 2 with independent parameter 
n. 

Computer program HEJET is used to perform the analysis described above. 

The results obtained from HEJET are provided in section 3.1. 

Considering the fact that similarity criteria is not satisfied under 
certain flow conditions, additional analysis is done with turbulent boundary 


u - U m f l<"> ’ T - T a = < T m ' T a> f 2<1> 


3 


n = r/(r&)u 
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layer code GENMIX (6). This program solves parabolic differential equations 
evolving from mass fraction, momentum and energy balance equations for single 
or multi -component system. Prandlt mixing length model is used for turbulence 
modeling with single equation. Analysis using GENMIX has also been presented 

X 

in section 3.2 and compared to HEJET for the particular cases analyzed here. 



3.1 Analysis using HEJET program 


Two sample cases have been chosen for calculation with HEJET program. The 
upstream conditions corresponding to these two cases are shown in Table 3.1. 
fiass flow rate of 0.8 kg/sec at a pressure of 1 bar is the common criteria for 
the two cases where the jet exit temperature is allowed to be 12 K and 120 K, 
respectively. 

Table 3.1 Helium Jet Properties used by HEJET Program 




Jet 

Jet 


Jet Velocity (m/sec) 

Pressure (bar) 

Temperature (k) 

Case 1 

10.82 

1 

12 

Case 2 

108.24 

1 

120 



/, 3 * 


Axial temperature profile for the two cases are illustrated in figure 3.1. 

The asymptotic value of temperature approaches the surrounding air temperature 
of 299k. While the temperature for case 2 reaches the asymptotic value before 
9 meters, corresponding temperature of case 1 is lower at the same axial 
distance. 


The potential core is characterized by the axial distance where/jet v e 1 o- - 
- city is equal to discharge velo city and temperature is equal to jet tem- 
perature. The potential core calculated for cases 1 and 2 are 1.38 and 0.94 
meters, respectively. Radial temperature distribution for cases 1 and 2 are 
shown in figures 3.2a and 3.2b respectively. Both cases show sharp gradient 
in temperature profiles which equilibrate to surrounding temperature of 299k 
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I 


I 


I 


at a normalized radial distance of 1. Radial profiles of temperature at 
several axial locations are shown in these figures. 

The jet centerline temperature and veloci ty profiles contrib ute the maxi- 

i " 'l 

Tnum radial magnit ude. .The radial distribution profiles begin^a normalized 

has been achieved from 


the relationship. 


Since (r&) u represents the radius of U max /2, which is linear with z, the 
profiles are regular. The initial assumption of Gaussian distribution of tem- 
perature and velocity reappear through the radial profiles. 

Axial velocity profile for cases 1 and 2 are shown in Figure 3.3. Very 
sharp decay of the velocity profile is evident from this figure for case 2. 

The asymptotic value of the velocity reaches the surrounding air velocity 
(which is stagnant in the present case). Radial velocity distribution is 
illustrated in figures 3.4a and 3.4b. 



distance 




The normalization of radial distance 


Radius of T_, /2 signifies the spreading rate of the expanding jet 
temperature profile. Figure 3.5 shows the radial T max /2 width of the jet in 
terms of axial distance. The linearity in the profile is the criteria for 
which self similarity of the jet is assumed. Approximate cor rejjuti-on--between 
the radial distance (r£) and axial position Z is such^fhat (rJ4) = .105 Z. 

_ — ^ ^ ^ ^ JL. - — 

Figure 3.6 illustrates radial density distribution function of case 1. The 
asymptotic value reached by the distributions is the surrounding air density 
at 299k. The actual value of helium gas density at jet discharge is greater 
than the surrounding air density (p^/p,,,. = 3.34) in case 1. Hence for the 

J c L all 
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mass flow rate of 0.8 kg/sec the jet velocity is low, which is equal to 10.824 
m/sec. On the other hand, for case 2, P: ai ./p^ tr . = 0.334 and velocity at 

jet a i r 

discharge equals to 108.24 m/sec. The asymptotic behavior of density for the 
two cases are therefore different as observed in figures 3.6a and 3.6b. 


% 

However, radial concentration distributions are asymptotically decreasing 
towards 0 at the jet periphery. Radial concentration distribution for cases 1 
& 2 are illustrated in'figures 3.7a and 3.7b respectively. 


Numerical values of the radial profiles are given in Table A.l as a 
function of axial distance and normalized radial distance. The normalization 
in any particular axial location is based on 3*(r5$) u , such that appropriate 
emphasis is given to quantities at any Z location. 
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3.2 Analysis using GENMIX program 


Gaseous helium jet dispersion to atmosphere has also been analyzed using 
the boundary layer code called GENMIX{6). Heated turbulent jet with com- 
bustion is typically analyzed using this program. As an effort to estimate 
the predictability of HEJET program; an analysis of cryogenic helium jet is 
presented here. Axial temperature and velocity calculations from the two 
programs are then compared. 


Input to GENMIX program was adjusted such that, 

XLAST = 9.5 meters; XOUT = 0.0; XEND = 0.0 

R b =0.0 ; R c = 0.0; = 0.0762 (jet radius) 

Calculation has been performed for helium gas using the following 
constants (5). 


R = 8314 J/kmol.k 

R„ = R/4.O03 ; C = - R ; C DUq = 5196.25 — 

He P 2 ’ p He Kg-K 

The same two reference cases as shown in section 3.1 for HEJET analysis is 

performed with GENMIX program, as shown in Table 3.2. 
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Table 3.2 •Helium jet conditions used by GENMIX program 



U jet (m/s) 

Tjettfc) 

Pressure (bar) 

Case 1 

10.82 

12 

1 

Case 2 

108.24 

120 

1 


Comparison of axial mean temperature profiles obtained for GENMIX and HEJET 
programs for Cases 1 and 2 are shown in figure 3.8. Axial temperature profi- 
les for case 1 are not in close agreement between the two analyses. Maximum 
deviation between the two calculations is approximately 40 K. GENMIX computer 
program is oriented towards parabolic solution which are dependent upon the 
initial conditions. The usual application of this code has been made to high 
temperature gaseous jets (6), where asymptotically decreasing temperature pro- 
files are traditionally analyzed. The case considered here has not been ana- 
lyzed jet by. GENMIX for validation against experiments. Hence, accuracy of 
either program in cryogenic temperature predictions remain to be performed. 

As observed from figure 3.8, the asymptotic temperature calculated by GENMIX 
at 10 meters appear to 80 K below the surrounding temperature. The possibi- 
lity of underprediction by GENMIX can only be confirmed by comparison with 
appropriate data. Radial temperature distribution from GENMIX analysis 
obtained for cases 1 and 2 are depicted in figures 3.9 and 3.10 respectively. 
Good agreement in axial velocity profiles of cases 1 and 2 are observed from 
figure 3.11. Corresponding radial velocity profiles are illustrated in 
figures 3.12a and 3.12b. . 
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The above analysis using two different computer programs provide a strong 
validity to calculation procedures. Prediction of helium jet dispersion to 
atmosphere can hence be continued using unsteady upstream conditions as input 
^to HEJET following the experimental conditions. Direct comparison of predic- 
tions by HEJET to actual available measurements is being proposed here in 
order to improve the working models for higher accuracy. 
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APPENDIX A 


User's guide for jet analysis program HEJET . 

This program calculates the centerline profiles and radial distribution 
profiles of 

1. temperature (or density, concentration) 

2. velocity 

of axisymmetric jets released to atmosphere. The solution scheme is based 
upon integral method and iterative approach as illustrated by References 1 & 

3. 

a) File management 

The main program will open three files for data input and output. These files 

Input data file 

Axial profiles of temperature, velocity, concentration 
Radial profiles of temperature, velocity, concentration 
and density 

Features of the above files are given below in detail: 

DATA. IN (Input data file) 

(A) ICASE: 0 = Indicates steady state calculation 
1 = Unsteady state calculation 

Default ICASE = 0, which implies steady state upstream 
condition is fixed by one set of input in card group (E). 


DATA-IN 
JET. OUT -#■ 

RADIAL. OUT 
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(B) DIA; ZTOT ; TTOT ; 


Diameter; total downstream distance; total time 


Diameter corresponds to pipe diameter at jet release 
(meters). Total downstream distance ZTOT is the 
distance through which calculation should proceed (m). 
Total time of unsteady upstream condition is to be spe- 
cified by TTOT. Default is TTOT = 0, which indicates 
steady state calculation. 


(C) Calculational Constants 


rr , P, , K e Ratio of 


, ir , entrainment constant 


Az , At , Ar : 

e c , KOUNT 
RMID, CMAX : 


Axial increment, temporal increment, radial 
increment 

convergence criteria, limit of iteration 
initial values of (rfc)u and C . 


(D) Surrounding air data 

RATM, TATM : Air density. Temperature. 


(E) Jet properties (upstream condition) 

This input should follow the sequence of ICASE, such that for i = 1, ITMAX 
the following data are given where ITMAX = ICASE + 1 = TTOT/DT 


FLOW (i), P(i), T(i): Mass flow rate (kg/sec). Pressure (bar). 

Temperature K 

End of file DATA. IN 
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JET. OUT 


(Output of axial profiles of temperature, velocity, 
concentration and density) 

(A) Characteristic numbers 
* Reynolds number 

Froude number 
Grashof number 

Pcore : Potential core length (m) 

(B) Axial output: Parameters as a function of axial distance at each time 

step [subscript max corresponds to centerline values]. 

OUTPUT Z :: R*. Re*. C max . Rho max , T(nax , U max 

Z : axial position 

R u : radius of U m /2 

n max 

RCw : radius of T „ v /2 

n max 

C max : Centerline concentration profile 

Rh°m ax : Centerline density profile 

T max : Centerline Temperature profile 

U m ax : Centerline Velocity profile 

End of file JET. OUT 

RADIAL. OUT -*• Output of radial profiles of temperature, velocity and 
concentration at each time step. 

For 0 < t < TTOT 

C ^ z < ZTOT 

Cs r n s 1 
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T (Z, r ) Temperature 

U (Z, r ) Velocity 

C (Z, r ) Concentration 

n' 

R (Z, r ) Density 

< 11 

here Z = axial distance 

r = normalized radial distance 
n 

■'» r / (3 * (R&)) 

End of file RADIAL. OUT 

Complete listing computer program HEJET is given in Appendix B. Output of 
sample cases I (discussed in Chapter 3) has also been provided in Appendix B. 

> 

i 
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FUNCl 


HEJET 


_ DATA. IN 
_ JET. OUT 
_ RADIAL. OUT 

-AINTG1 


RENUMB 

FNUMB 


MIDPNT 


-QROMO (3) — 


POLINT 


FUNC2 

FUNC3 


Table A.l Structure of subroutines and I/O files of Helium jet 
dispersion code HEJET 
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velocity 


Output of HEJET p r ogram used to calculate the 
axial and radial distribution of temperature and 

of a gaseous Helium jet dispersed into atmosphere. 


CASE 1 


OUTPUT. FILE : JET. CUT 


************* jet input axrornoNS ************* 

Time step Flew rate(kg/sec) Pressure (bar) 

Temperature Viscosity 

1 0.800000,01 1.00000000 12.00000000 
2 . 23000006E-06 

**************** COTPOT OF HEJET ****************** 
******************* Dimensionless Groups********* 
renolds nunber= 2 . 99772775E+06 fraud nunto= 112.32895700 

grashof number= 3.74713672E-05 pcore = 1.37915623 

Jet Density Jet Velocity jet Temperature 
4.05232 10.82452 12.00000 


*********** Axial distributions 


z 

Rl/2 

Rcl/2 

Unax 

1.000000 

0.066777 

0.078129 

10.824515 

2.000000 

0.164225 

0.192143 

8.539839 

3.000000 

0.258746 

0.302732 

5.911030 

4.000000 

0.352793 

0.412768 

4.508695 

5.000000 

0.446674 

0.522608 

3.641720 

6.000000 

0.540477 

0.632358 

3.053618 

7.000000 

0.634236 

0.742056 

2.628735 

8.000000 

0.727972 

0.851728 

2.307508 

9.000000 

0.821691 

0.961379 

2.056164 

10.000000 

0.915398 

1.071016 

1.854156 

Number of iterations 

= 40 


************** 


Qnax 

Rhcmax 

Tmax 

1.225565 

4.052319 

12.000000 

0.695905 

2.381118 

99.275238 

0.478125 

1.838207 

161.778152 

0.363496 

1.641240 

194.676529 

0.293055 

1.539845 

214.893234 

0.245436 

1.478115 

228.559860 

0.211110 

1.436600 

238.411407 

0.185199 

1.406775 

245.847824 

0.164949 

1.384314 

251.659683 

0.148688 

1.366791 

256.326538 






OOTEUr FILE: RADIAL. CUT 


Radial Distribution of Temperature 


12.000 12.000 

65.478 

129.222 

186.317 

230.725 

261.235 

279.930 

290.209 295.301 

*100 ..769 116.477 

145.577 

181.270 

216.527 

246.258 

268.208 

282.589 

291.015 295.453 

162.193 171.319 

190.215 

214.388 

238.921 

260.056 

275.955 

286.551 

292.860 296.236 

194.846 201.184 
294.086 296.774 

215.137 

233.362 

252.101 

268.409 

280.784 

289.098 

214.979 219.805 
294.919 297.145 

230.856 

245.472 

260.615 

273.871 

283.982 

290.807 

228.609 232.496 
295.516 297.412 

241.640 

253.836 

266.536 

277.698 

286.239 

292.021 

238.442 241.691 
295.963 297.613 

249.489 

259.951 

270.885 

280.521 

287.912 

292.926 

245.868 248.656 
296.309 297.770 

255.453 

264.613 

274.211 

282.686 

289.199 

293.625 

251.674 254.115 

260.138 

268.283 

276.835 

284.400 

290.220 

294.180 

296.584 297.895 

256.337 258.507 

263.914 

271.246 

278.959 

285.788 

291.049 

294.632 

296.809 297.997 




Radial Distribution of Velocity 



10.172 

0.069 

8.438 

0.021 

6.179 

3.994 

2.279 

1.148 

0.510 

0.200 

8.453 

0.104 

7.549 

0.034 

5.952 

4.142 

2.545 

1.380 

0.661 

0.279 

5.887 

0.084 

5.356 

0.028 

4.301 

3.049 

1.908 

1.054 

0.514 

0.221 

4.499 

0.069 

4.128 

0.023 

3.344 

2.391 

1.509 

0.841 

0.414 

0.180 

3.637 

0.058 

3.354 

0.020 

2.730 

1.962 

1.244 

0.697 

0.344 

0.150 

3.051 

0.050 

2.822 

0.017 

2.305 

1.662 

1.058 

0.594 

0.295 

0.129 

2.627 

0.044 

2.436 

0.015 

1.994 

1.441 

0.919 

0.517 

0.257 

0.113 

2.306 

0.039 

2.142 

0.013 

1.756 

1.271 

0.812 

0.458 

0.228 

0.100 

2.055 

0.035 

1.912 

0.012 

1.569 

1.137 

0.728 

0.411 

0.205 

0.090 

1.854 

0.032 

1.726 

0.011 

1.418 

1.029 

0.659 

0.372 

0.186 

0.082 


Radial Distribution of Concentration 
0.956 0.834 0.664 0.483 0.320 0.194 0.107 0.054 

0.025 0.011 



0.564 

0.023 

0.519 

0.010 

0.436 

0.335 

0.234 

0.150 

0.088 

0.047 

0.389 

0.017 

0.363 

0.008 

0.309 

0.241 

0.171 

0.111 

0.066 

0.035 

0.296 

0.014 

0.278 

0.006 

0.238 

0.187 

0.133 

0.087 

0.052 

0.028 

0.239 

,,’0.012 

0.225 

0.005 

0.194 

0.152 

0.109 

0.071 

0.043 

0.023 

0.200 

0.010 

0.189 

0.005 

0.163 

0.128 

0.092 

0.061 

0.036 

0.020 

0.172 

0.009 

0.163 

0.004 

0.141 

0.111 

0.080 

0.053 

0.032 

0.017 

0.151 

0.143 

0.124 

0.098 

0.070 

0.046 

0.028 

0.015 

0.008 

0.003 







0.135 

0.007 

0.128 

0.003 

0.110 

0.087 

0.063 

0.042 

0.025 

0.014 

0.121 

0.006 

0.115 

0.003 

0.100 

0.079 

0.057 

0.038 

0.023 

0.012 




Radial * Distribution of Density 



4.052 

1.252 

4.052 

1.236 

2.834 

2.086 

1.687 

1.469 

1.349 

1.285 

2.364 

1.249 

2.202 

1.236 

1.954 

1.716 

1.532 

1.405 

1.324 

1.276 

1.835 

1.244 

1.776 

1.233 

1.665 

1.542 

1.435 

1.353 

1.298 

1.263 

1.640 

1.240 

1.607 

1.232 

1.539 

1.458 

1.383 

1.323 

1.282 

1.255 

1.539 

1.237 

1.517 

1.231 

1.468 

1.408 

1.351 

1.305 

1.271 

1.250 

1.478 

1.235 

1.461 

1.230 

1.423 

1.376 

1.330 

1.292 

1.264 

1.246 

1.436 

1.234 

1.423 

1.229 

1.393 

1.353 

1.315 

1.283 

1.259 

1.243 

1.407 

1.233 

1.396 

1.229 

1.370 

1.337 

1.304 

1.276 

1.255 

1.241 

1.384 

1.232 

.1.375 

1.228 

1.353 

1.324 

1.295 

1.270 

1.252 

1.240 

1.367 

1.232 

1.359 

1.228 

1.339 

1.314 

1.288 

1.266 

1.249 

1.238 


CASE 2 


I. OUTPUT FITE : JET. OUT 


************* jet INPUT OCNDITICNS ************* 







Time step Flow rate(kq/sec) Pressure (bar) 

Temperature Viscosity 

1 0.80000001 1.00000000 120.00000000 

1.14000002E-04 

**************** OUTPUT OF HEJET ****************** 

* 

************** Dinensionless nr mip^ ********* 

renolds nuniber= 58639.76170000 froud numb= 3874.14185000 
grashof number^ 6.60668090E-02 pcore = 0.94162774 

Jet Density Jet Velocity jet Temperature 
0.40523 108.24515 120.00000 

*********** Axial di stri b utions ************** 


z 

Umax 

Rl/2 

Rcl/2 

Qnax 

Rhomax 

Tmax 

0.500000 

108.245148 

0.061999 

0.072539 

1.112296 

0.376908 

120.000000 

1.500000 

43.170197 

0.157886 

0.184727 

0.339863 

0.725913 

238.164536 

2.500000 

25.150795 

0.252165 

0.295033 

0.199064 

0.873316 

263.367493 

3.500000 

17.716330 

0.346118 

0.404958 

0.140582 

0.953759 

273.835815 

4.500000 

13.667797 

0.439950 

0.514742 

0.108620 

1.004317 

279.557037 

5.500000 

11.123327 

0.533726 

0.624459 

0.088486 

1.039012 

283.160919 

6.500000 

9.376678 

0.627466 

0.734136 

0.074644 

1.064290 

285.638672 

7.500000 

8.103730 

0.721187 

0.843789 

0.064545 

1.083523 

287.446503 

8.500000 

7.134884 

0.814895 

0.953427 

0.056851 

1.098647 

288.823669 

9.500000 0.908593 

6.372854 

Number of iterations 

1.063054 
= 28 

0.050796 

1.110852 

289.907593 


OUTPUT FILE : RADIAL. CUT 


Radial Distribution of Temperature 


120.000 135.374 
294.329 297.048 
238.657 243.521 
296.598 297.935 
263.481 265.871 
297.414 298.286 
273.878 275.415 
297.818 298.465 


169.599 205.579 
252.435 263.321 
270.791 277.073 
278.786 283.184 


237.429 261.954 
274.043 283.063 
283.440 288.920 
287.703 291.634 


278.652 288.797 
289.710 294.056 
293.039 295.782 
294.615 296.617 


279.577 280.697 
298.058 298.572 
283.172 284.048 
298.217 298.643 
285.646 286.363 
298.331 298.695 
287.451 288.058 
•298.415 298.733 
288.827 289.352 
298.481 298.763 
289.910 290.373 
298.533 298.786 


283.254 286.634 
286.106 288.849 
288.084 290.392 
289.536 291.527 
290.647 292.398 
291.525 293.087 


290.134 293.197 
291.705 294.213 
292.803 294.927 
293.613 295.455 
294.237 295.862 
294.731 296.186 


295.533 297.109 
296.133 297.432 
296.556 297.662 
296.871 297.832 
297.113 297.964 
297.306 298.070 


Radial Distribution of Velocity 


100.714 

0.636 

82.747 

0.193 

60.012 

38.419 

21.711 

10.831 

4.769 

1.854 

42.693 

0.517 

38.053 

0.170 

29.940 

. 20.794 

12.749 

6.899 

3.296 

1.390 

25.041 

0.355 

22.764 

0.119 

18.267 

12.939 

8.091 

4.466 

2.176 

0.936 

17.675 

0.269 

16.212 

0.091 

13.127 

9.382 

5.919 

3.296 

1.620 

0.703 

13.648 

0.217 

12.583 

0.074 

10.240 

7.356 

4.665 

2.611 

1.290 

0.563 

11.113 

0.181 

10.279 

0.062 

8.393 

6.050 

3.849 

2.162 

1.072 

0.469 

9.370 

0.156 

8.688 

0.053 

7.110 

5.137 

3.276 

1.844 

0.916 

0.402 

8.099 

0.136 

7.522 

0.047 

6.167 

4.463 

2.851 

1.608 

0.800 

0.352 

7.132 

0.121 

6.633 

0.042 

5.445 

3.946 

2.524 

1.425 

0.710 

0.313 

6.371 

0.109 

5.931 

0.038 

4.874 

3.536 

2.264 

1.280 

0.639 

0.281 




Radial Distribution of Concentration 


0.949 

0.023 

0.822 

0.010 

0.650 

0.469 

0.309 

0.186 

0.102 

0.051 

0.303 

0.012 

0.279 

0.005 

0.234 

0.179 

0.125 

0.080 

0.047 

0.025 

0.178 

0.008 

0.166 

0.004 

0.142 

0.110 

0.078 

0.051 

0.030 

0.016 

0.126 

0.006 

0.118 

0.003 

0.102 

0.079 

0.057 

0.037 

0.022 

0.012 

0.098 

0.005 

0.092 

0.002 

0.079 

0.062 

0.045 

0.029 

0.017 

0.009 

0.079 

0.004 

0.075 

0.002 

0.065 

0.051 

0.037 

0.024 

0.014 

0.008 
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0.067 

0.003 

0.063 

0.002 

0.055 

0.043 

0.031 

0.020 

0.012 

0.007 

0.058 

0.003 

0.055 

0.001 

0.048 

0.038 

0.027 

0.018 

0.011 

0.006 

0.051 

0.003 

0.048 

0.001 

0.042 

0.033 

0.024 

0.016 

0.009 

0.005 

0.046 

* 6.002 

0.043 

0.001 

0.038 

0.030 

0.021 

0.014 

0.009 

0.005 




Radial Distribution of Density 



0.391 

1.164 

0.430 

1.199 

0.497 

0.596 

0.722 

0.863 

0.996 

1.098 

0.728 

1.193 

0.753 

1.210 

0.803 

0.873 

0.956 

1.038 

1.109 

1.160 

0.874 

1.203 

0.891 

1.215 

0.929 

0.982 

1.042 

1.100 

1.148 

1.182 

0.954 

1.209 

0.967 

1.218 

0.997 

• 1.039 

1.086 

1.131 

1.167 

1.193 

1.005 

1.212 

1.015 

1.219 

1.040 

1.075 

1.113 

1.150 

1.179 

1.199 

1.039 

1.214 

1.048 

1.220 

1.069 

1.099 

1.132 

1.162 

1.187 

1.204 

1.064 

1.216 

1.072 

1.221 

1.090 

1.116 

1.145 

1.171 

1.192 

1.207 

1.084 

1.217 

1.090 

1.221 

1.107 

1.130 

1.155 

1.178 

1.196 

1.209 

1.099 

1.218 

1.105 

1.222 

1.119 

1.140 

1.162 

1.183 

1.199 

1.211 

1.111 

1.116 

1.130 

1.148 

1.169 

1.187 

1.202 

1.212 


1.219 1.222 


PROGRAM HEJET 


c 

c************************************** ************************** 
c this program calculates the centerline decay profiles 

c and radial distribution profiles 

c 1. concentration { density , temperature) 

c 2. velocity 

,c of axisymmetric jets by integral method 

c Reference : Chen and Rodi ( Pargamcan press) 

c William Pitts ( NBS report ) 

c*** ************************************************************ 

c 

parameter ( jmax = 10 , kmax= 10 ) 
c 

dimension flow( 20 ), p(20) , t(20) , uo( 20) , ro(20) , 

vis (20) 
c 

dimension ru2(kmax, 20 ), uc(kmax, 20 ), oc(kmax, 20) 

1 , tc(kmax, 20 ), rc(kmax, 20 ) 

dimension cr(jmax, kmax, 20 ) , ur(jmax, kmax, 20 ) 

1 , tr(jmax, kmax, 20 ) , dr(jmax, kmax, 20) 

2 , cm(jmax, kmax, 20 ) 

c 

data roans , peons , agrav / 2.077 , 1.01e2 , 9.81 / 

c*************************************************************** 
open (unit = 8 , file= / c:\nunirec\data.in / ,status== , old / ) 
epen (unit = 6 , 

filep= / c: VJ0MREC%j et . cut ' , status= / uriknown' ) 
open (unit = 7 , 

file='c: \NUMREC\radial . cut' , status= / unknown / ) 

c* ************************************************************** 

C DATA ETTTRY 

c Specify steady state or transient clculation : 

c 

read(8,*) icase 
c 

C A) Read geometry 

c 

read(8,*) dia, zin, ztot, ttot 
C 

c B) Read calculation constants 

c 

read(8,*) rr, pi, ake 

read(8,*) dz, dt, delr 

read(8,*) epsc, kount 

read(8,*) rmid, anax 

c 

c C) Read atmospheric data 

c 

read (8,*) ratm, tatm 
c 



+> . 


c D) Read jet properties 
c 

itmax = 1 

c iff icase. eq. 0) itmax = 1 

if ( icase. eq. 1) itmax = int( ttot/dt ) 
c if(icase.eq. 1) itmax = 2. 

c 

write(6,*) ************** JET INHJT CONDITIONS 
************** 

write(6,*) 'Time step Flow rate (kg/sec) Pressure (bar) 

1 Temperature Viscosity ' 
do 10 it = 1, itmax 

read(8,*) flow( it ), p( it ) ,t ( it ) , vis( it ) 
write(6,*) it, flcw( it) , p(it) , t(it) , vis(it) 

10 cxoTtinue 

write(6,*) ***************** OUTFOT OF HEJET 
****************** • 

C 

c Conversion to density and velocity at the inlet 

c 

area = pi * dia**2. / 4. 

do 20 it = 1, itmax 

po = p( it ) * paons 

to = t( it ) 

vsp = icons * to / po 

dens = 1./ vsp 

mo ( it'.) - dens 

uzero = flew (it) / ( area * ro(it) ) 

uo ( it ) = uzero 
20 continue 

c*************************************************************** 

c 

c Begin solution 

c Calculate the dimensionless groups 

c 

do 300 it = 1, itmax 
c 

write (6 ,*) ******************** Dimensionless 
Groups********* # 

dens = ro (it) 
uzero = uo (it) 
vise = vis (it) 

reno = mumb( uzero , dens, dia, vise ) 

froud = fnurab( uzero , dens, dia, ratm, agrav) 

grash = froud / reno 
xc = 2.13 * reno **.097 

pcore = dia * xc 

write(6,*) 'renolds number= ' ,reno, ' froud numb= ' , 

fraud 

write (6,*) 'grashof number= ' , grash, ' pcore = pcore 
c 

rjet = ro( it ) 

fac = ro( it ) / ratm 



z = zin 

ror = dia/ 2.0 

reps = nor * ( ro(it)'/ ratra )** 0.5 

writers,*) 7 Jet Density ' , 7 Jet Velocity 7 jet 
Temperature 7 

write(6 f 30) rjet, uo(it) , t(it) 
c 

x'30 format(3(4x,fl0.5) ) 
kk = 1 

c begin iterative solution for each z location 
c 

write(6,*) 7 *********** Axial distributions 
************** > 

write(6,*) 7 z Rl/2 Rcl/2 Cmax 

Rhcanax 

1 Tmax Umax 7 


c 

do 200 k = 1, kmax 
1 cmaxo = cmax 
c . 

c calculate the integrals airrtl aint2 & aint3 
call aintgl ( rmid, anax, rjet, ratm,rr 
1 , aintl, aint2, aint3 

) 


c 

c 

c 


c 

c 

c 

c 


c 

rjet 


calculate cmax 


cofl = ( pi**0.5) * reps / ( ake * z ) 
cmax = ( aint3 / aint2 ) * oofl 

dele = cmax - cmaxo 


if ( abs(delc) -lt.epsc*cmax) then 


Centerline decay profiles 
oc(k, it) = cmax 


cof2 

ru2(k,it) 

rc2 

cof3 

uc(k,it) 

rc(k,it) 

tc(k,it) 


= ( ate * z ) / ( 2. * pi )**0.5 
= (cof2 * aintl**0.5 ) / airrt3 
= rr * ru2 (k, it) 

= ( pi**0.5) * reps / ( ate*z) 

= uo(it) * cof3 * ednt3 / aintl 
= ratm /( 1.+ cmax* ( l./fac - 1.)) 
= tatm + cmax*( t(it) - tatm ) 


if ( tc(k,it) .It. t(it) ) tc(k,it) = t(it) 

if( rjet .gt. ratm. and. rc(k,it) .gt. rjet) rc(k,it) = 


if ( uc(k,it) .gt. vizero ) uc(k,it) = uzero 


c 

write(6,90 ) z, ru2(k,it), rc2, cc(k,it), rc(k,it), 

tc(k,it) 

1 , uc(k,it) 

90 format ( 7(lx,fl0.6) ) 


c 


else 

if (kk.eq.kount) go to 1000 
kk = kk + 1 
go to 1 

endif 

c 

c***************************************************** 

<£• calculate the radial profiles for variables 
c ur = velocity, cr= concentration 

c dr = density , tr= temperature 

c **************************************************** 

C 

r = .02 
c 

delr = (ru2(k,it) * 3. ) / float (j max) 

do 100 j =1, jmax 
c 

if( k.eq.l. and.j.eq.l ) cznorm = cc( k, it ) 
ratio = r/ ru2 (k,it) 

power = HD. 693 * ratio* *2. 

ur (j,k,it) = uc(k,it) * exp (power) 

cr (j,k,it) = cc(k,it) * exp(pcwer/rr**2.) 

cm(j,k,it) = cr^jfkfit) / cznorm 

dr (j,k,it) = ratm / ( 1 + cr(j,k,it)*( rataVrofit) -1 

)) 

tr (j,k,it) = tatm + cr(j,k,it)*( t(it) - tatm ) 

r - r + delr 

c 

if ( ur(j,k,it) . gt. uzero ) ur(j,k,it) = uzero 
if( tr(j,k,it) . It. t(it) ) tr(j,k,it) == t(it) 
if( rjet .gt. ratm. and .dr(j,k,it) .gt.rjet) 
dr(j,k,it)=rjet 
c 

100 continue 
c 

c if( k.eq. 1) z = zin 

z = z + dz 
c 

200 continue 
c 

300 continue 
c 

c write the radial profiles 

c 

wrrte(7,*)' Temperature ' , ' Velocity 
Concentration ' 

1 , ' Density ' 

c 

do 92 jj = 1, itmax 
format (10(lx,f7.3)) 

write(7,91) (( tr (ii,kz,jj) , ii= 1, jmax) ,kz= 1, kmax) 

write(7,91) (( ur (ii,kz,jj) , ii= l, jmax),kz=l, kmax) 

write(7,91) (( cxn(ii,kz, jj) , ii= 1, jmax) ,kz= 1, kmax) 


91 



write(7,91) ( ( dr (ii,kz,jj) , ii= 1, jmax) ,kz= 1, kmax) 

92 continue 

c 

1000 write (6, 1001) kk 

1001 format (2x, ' Number of iterations = ',15) 

close(8) 

x’ close(6) ' - 

close(7) 
end 

function rnumb(u, r, d, v ) 

muxnb = u * r* d/ v 

end 

Q^-k'k-k'k-kick-kicIck-k-k'kifk-k'k-k-k-klrklck'k-kitirk-k-k-k-k'k-k-k'k-k-k'k-k-ick'k'klckick-k-k-k'k-k'k-k 

function fnumb( u, r, d, rat, g ) 
artum = u **2.0 

deno = g * d * abs( rat - r ) / r 
fnumb = anum / deno 
end , 

c 

SUBROUTINE MDFNT{ A, B, HYPE, S, N ) 

c 

c Ihis routine computes tne n'th stage of refinement of 

an extended 

c midpoint rule. FUNC is ir?xrt as the name of the 

function 

c to be integrated between limits A and B . When called 

by 

c N = 1 , the noutien returns as S the crudest estimate 

of 

c estimate of int f(x)dx. As N increases the accuracy 

c increases by ( 2/3 )* 3**(N-1) additional interior 

points, 
c 

if ( itype.eq.l) then 
C 

if (n . eq. 1 )then 

s = (b-a) * funcl( 0.5* (a + b)) 
ifc= 1 
else 

tnm = it 

del = ( b-a ) / (3. * tnm ) 
ddel = del + del 
x = a + 0.5 * del 
sum = 0.0 
do 11 j= l,it 

sum = sum + funcl(x) 
x = x + ddel 
sum = sum + funcl(x) 
x = x + del 
11 continue 
C 

s = ( s+ ( b-a ) * sum/ trim) / 3. 



it = 3 * it 

endif 
c 

elseif( itype. eq. .2 ) THEN 
c 

if (n . eq. 1 )then 

s = (b-a) * func2( 0.5* (a + b) ) 

V it= 1 

else 

tnm = it 

del = ( b-a ) / (3. * tnm ) 
ddel = del + del 
x = a +. 0.5 * del 
sum =0.0 
do 22 j= l,it 

sum = sum + func2(x) 
x = x + ddel 
sum = sum + func2 (x) 
x = x + del 
22 continue t 

C 

s = ( s+ ( b-a ) * sun/tnm) / 3. 

it = 3 * it 
endif 
C 

else 

C 

if (n . eq. 1 )then 

s = (b-a) * func3( 0.5* (a + b) ) 
it= 1 
else 

tnm = it 

del = ( b-a ) / (3. * tnm ) 
ddel = del +.del 
x = a + 0.5 * del 
sum = 0.0 
do 33 j= l,it 

sum = sum + func3 (x) 
x = x + ddel 
sum = sum + func3(x) 
x = x + del. 

33 continue 
C 

s = ( s+ ( b-a ) * sum/ tnm) / 3. 

it = 3 * it 

endif 
C 

endif 

C 

return 

end 

C***************************************************** 

SUBROUTINE FOUNT ( XA, YA, N,X, Y,DY) 


C 



] 


4 


I 


c 

C 

C 

c 

c 

C 

< ' 
C 


11 


c 

c 


12 


13 


Given arrays XA and YA , each of length N, 
and given a value of X , this routine 
returns a value Y and an error estimate 
DY. If P(x) is a polynomial of degree 
N-l , then it returns Y= P(X) 

parameter (rnnax = 10) 

dimension xa(n) , ya(n) , c(nmax) , d(nmax) 
ns =1 

dif = abs( x - xa(l) ) 
do 11 i = 1, n 

dift - abs( x- xa(i) ) 
if ( dift. It. dif) then 
ns = i 
dif = dift 
endif 

c(i) = ya(i) 
d(i) = ya(i) 
continue 

= ya( ns) 

= ns - 1 

m= 1, n-l 

i =1, n-ffl 
ho = xa(i) - x 

l?) = xa(i+m) -x 

w = c( i+l) - d(i) 

den = ho - hp 

if ( den. eq. 0.0) pause 
den = w / den 
d(i)= hp * den 
c(i)= ho * den 
cxxitinue 

if ( 2*ns. It. n-m ) then 


dy 

= c(ns + l) 

else 



= d(rss) 

ns 

= ns -1 

_ endif 


Y 

- Y + DY 

cxxxtinue 

return 

end 



Y 

ns 

do 13 
do 12 


^k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k’k'k'k'k'k'k'k’k'kit 


c 

c 

c 

c 

c 


SUBROUTINE QRCMD( A, B, HYPE, SS, DSS ) 

Romberg integration on an open interval. 

Returns by SS the integral of the function FUNC 
f ro m A to B, using Subroutine MUJFNT 


parameter (eps = 1. 

k = 7 


e-5 


jmax= 15, jmaxp=jmax + 1 
km = k- 1 ) 


1 


9 


/ 



) 


> 


c 

dimension s( jmaxp), h( jmaxp) 
c 
C 

h(l) =1. 
do 11 j = 1, jmax 

call midpnt( a, b, itype, sp, j ) 

V s(j) = sp 

if ( j. ge. k) then 

call polint( h(j- km), s(j- km), k,0.,ss,dss) 

if (abs( dss).lt. eps* abs(ss) ) return 
endif 

s(j+l) = s(j) 
h(j+l) = h(j) / 9. 

11 continue 

write(*,*) j,sp,s(j-2),s(j-l),s(j) 

pause ' too many steps' 

end 

C********* ********************************************** 

c 

C 

SUBRCXJITNE aintgl(rmidl,onaxl, rjet,ratm,r 
1 , aintl, aint2, aint3 ) 

C 

C This program is only to test the integration 

c routine 

C 

ocmrcn/ccnd/rmid, cmax, rj, rat,rr 
nnid = rmidl 
cmax = anaxl 
rj = rjet 
rat - ratm 
rr = r 
c 

a = 0.0 
b = 2.0 
C 

itype = l 

call qncmo ( a, b, itype, ss, dss ) 

aintl = ss 
C 

itype = 2 

call qrano ( a,b, itype, ss, dss ) 
aint2 = ss 
C 

itype = 3 

call qrctno ( a, b, itype, ss, dss ) 

aint3 = ss 
C 

return 

HID 

C************************** *************************** 
FUNCTION FUNCl(XX) 



c 


c 

£ 


c 


cxsnrnon/cxaiil/ rmid ,anax, rjet, ratm, rr 

dcoef = cnax*(ratm / rjet -1.) 
dpow = -.693 /rr**2. . 
dpcwl = dpcw * (xx / rmid) **2 . 

denol = 1. + dcoef* exp( dpowl ) 

apcw = -.693* (xx/rmid) **2. 
anura = exp( apcw ) 

funcl = (( anum**2. / denol ) / rmid**2.) * xx 


C 


END 

Qlc-kicirklcirkltirk'klritirklticklilCklrklcIrklrkirtrklHrk-k-k-lrirtcirk-klcicklrkirk-k-k'k-kicirk'k-k 

JUNCTION FUNC2 (XX) 

OCMOT/OCM1/RMID , CMAX, RJET, RATM, RR 


dcoef = cmax*(^atm / rjet -1.) 
dpow = -.693 /rr**2. 
dpcwl = dpow * (xx / rmid) **2 . 


c 


denol — 1 . + dcoef* exp( dpowl ) 
c 

apcw = -.693* (xVrmid)**2. 
apcwl = apcw * ( 1. + 1./ rr**2. ) 
anum = exp( apowl ) 
c 

func2 = (( anum / denol ) / rmid**2.) * xx 
c 


END 

q*********************************************************** 

FUNCTION FUNC3 (XX) 

cem^/coo/emid , cmax, rjet, rajm, rr 


dcoef = cmax* (ratm / rjet -1.) 
dpow = -.693 /rr**2. 
dpcwl = dpcw * (xx / rmid) **2 . 


denol = 1. + dcoef* exp( dpcwl ) 

apcw = -.693* (xx/rmid) **2 . 
anum = exp( apcw ) 


c 


func3 = (( anum / denol ) / rmid**2 . ) * xx 
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Figure 3.3 
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Figure 3.4a Radial velocity profile of Helium jet at 12 K 
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Figure 3.9 Radial temperature distribution obtained from GENMIX 
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